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The exact treatment of Markovian open quantum systems, when based on numerical diagonalization of the 
Liouville super- operator or averaging over quantum trajectories, is severely limited by Hilbert space size. 
Perturbation theory, standard in the investigation of closed quantum systems, has remained much less 
developed for open quantum systems where a direct application to the Lindblad master equation is 
desirable. We present such a perturbative treatment which will be useful for an analytical understanding of 
open quantum systems and for numerical calculation of system observables which would otherwise be 
impractical. 

Understanding open quantum systems is of fundamental importance in many contexts of physics, ranging 
from traditional atomic and molecular physics 1 " 3 to recent studies of circuit QED 4 " 7 and optomechanical 
systems 8 " 10 . However, the simulation of open quantum systems is even more demanding than that of 
closed systems, and analytical and numerical approximation methods are less developed. This makes theoretical 
studies of large open systems particularly challenging and has spurred immense interest in the development of 
open-system quantum simulators 11 which could also be used to study dissipative phase transitions 12 " 14 . 
Markovian open quantum systems can be described by the Lindblad master equation 15 , 



p(0=Lp(0. 



(i) 



which governs the time evolution of the density matrix p in terms of the generalized Liouville super- operator L. 
Equation (1) allows for the investigation of system dynamics and of the asymptotic steady state by exact 
diagonalization of L or by simulating quantum trajectories 16 . However, the exponential growth of the Hilbert 
space dimension N with system size and the even more dramatic growth of the dimension N 2 of the associated 
space of linear operators severely limit the possibility of obtaining numerically exact solutions. Approaches 
alternative to the exact diagonalization of L include sophisticated numerics such as matrix product methods 17 " 19 
and self-consistent projection operator methods 20 . Although these techniques are suitable for larger systems, they 
are only easily applicable to one-dimensional or highly symmetric systems. 

For closed quantum systems of large size, perturbation theory is routinely employed as a useful approximation. 
In the past, specialized perturbative methods have also been applied to open quantum systems in specific contexts 
including full counting statistics 21 ' 22 , finite-time evolution 23 and critical behavior 1314 . One example of an open- 
system perturbative treatment is adiabatic elimination (generalized Schrieffer- Wolff formalism), which gives a 
simplified Liouville super-operator 24 " 27 whenever the spectrum separates into slow and fast degrees of freedom. 
Another method consists of a perturbative expansion of dissipative terms 23 ' 28 , which can yield the approximate 
finite- time evolution of systems with small dissipation rates. Moreover, a perturbative treatment for open systems 
can be formulated within the Keldysh Green's function framework based on bosonic or fermionic field 
propagators 1314,29 . 

Here, we present a canonical perturbative method that systematically determines the corrections to eigenstates 
and eigenvalues of the Liouville super- operator L Our treatment is not limited to a certain system or system type. 
Rather, it is applicable to a wide range of perturbations and different open quantum systems. Perturbation theory 
(PT) similar to the one presented here was previously studied 30 ' 31 but results were limited to the steady state and 
the issue of non-positivity of density matrices due to truncation was not addressed by the authors. 

The density-matrix PT we develop yields results both for the steady state as well as for all other eigenstates of L 
Further, we propose and derive a new PT based on the amplitude matrix. This amplitude-matrix PT guarantees a 
properly positive steady-state density matrix. Both kinds of PT can be applied to large open systems with lattice 



SCIENTIFIC REPORTS | 4:4887 | DOI: 1 0.1 038/srep04887 



1 



structure. Understanding the properties of such open systems is of 
immense interest and would otherwise be impractical due to the 
system size. The study of such systems is an important motivation 
for developing the PT presented here. Its application to open lattice 
systems is beyond the scope of the present paper and will be pub- 
lished elsewhere. 

Results 

Density-matrix PT. We propose a non-degenerate density-matrix 
PT based on the quantum master equation (1). The Liouville super- 
operator L, which serves as the generator of the quantum dynamical 
semigroup 15 , is generally not Hermitian, i.e. the adjoint super- 
operator is not equal to L The right and left eigenstates and 
of L associated with the eigenvalue k^eC are defined by 

lu^ = ^u^ L f w At = (^)*uy (2) 

Here, n is a non-negative integer labeling the eigenstates. Suitably 
normalized, the left and right eigenstates obey the bi-orthonormality 
relation (w^, u v ) = S^, where (x,y) = Tr [x f y] is the Hilbert-Schmidt 
inner product 32 for linear operators x and y. Together, the 
eigenvalues k^ and the associated right eigenstates contain all 
information of the steady state (labeled by \i = 0) and the 
dynamics of the system. 

In analogy to closed-system PT, the density-matrix PT is based on 
the series expansion of k^ and with respect to a small parameter a set 
by the perturbation. Starting point is, thus, the separation of L into two 
parts: an unperturbed super-operator L 0 and a perturbation a D_i, i.e. 

L = L 0 + aLi. (3) 

Like the original Liouville operator, L 0 must still be a proper gen- 
erator of the quantum dynamical semigroup. In addition, L 0 should 
be solvable in the sense that its unperturbed spectrum {k^\ 
u^} is known, or, at least a subset of interest. As appropriate for 
non-degenerate PT, we will assume that this part of the spectrum is 
non- degenerate. 

We next employ a general series expansion of eigenvalues and 
eigenstates in a, 

GO GO 

^=E^/, ty=E«^- (4) 

j=0 j=0 

Here, the index j counts orders of perturbation theory and k^ and 
ufjp are hence the j-th-order corrections to the eigenvalue and the 
right eigenstate. We determine recursion relations for k^ and by 
plugging equations (3) and (4) into equation (2) and examining the 
result order by order in a. For the j-th order in a, we obtain the 
general recursive expression 

(l 0 -4°)) M 0) = _ Ll „&-i) +A « (5) 

withA^^^*). 

So far, our treatment directly mirrors the well-known derivation of 
stationary PT for a closed system. Specifically, replace L 0 , li, and 
kfj, by the unperturbed Hamiltonian H 0 , the perturbation V, the 
eigenvectors and eigenenergies of H = H 0 + aV. Then, 
equation (5) takes exactly the form of the usual recursion equation 
in closed-system PT, namely 

(h 0 -EM)\^) = -v\^) + SU\ (6) 
with = Y^ k= i E f ] \^~ k) )- To obtain the 

recursive relation for 

the eigenenergy correction we multiply equation (6) with (x//^ \ 
from the left. This yields 



k=i 

Analogously, we take the inner product of equation (5) with the left 
eigenstate . This yields the recursion relation for k^ , 

A» = (^L 1 «0-'))-gAW(^),«0-*)>. (8) 

k=l 

Note that equation (7) is usually simplified further by demanding 
(x//^ | xj/^ ) = 0 for; # 0. We will instead keep the corresponding term 
(w^\ u^~ k ^) in equation (8) for reasons that will become clear 
momentarily. 

We next turn to the computation of the eigenstate corrections. The 
Hamiltonian of any closed system is Hermitian. Hence, H 0 provides a 
complete orthonormal eigenbasis {|<A^)}. As a result, in equa- 
tion (6) can be expanded in this eigenbasis. Solving equation (6) is 
then straightforward. By contrast, L is not Hermitian and may not 
even be diagonalizable. As a result, the expansion of uff in terms of 
will generally fail. We therefore adopt the different strategy of 
applying an appropriate generalized inverse to the singular super- 
operator (L 0 — k^). Several options exist for such a generalized 
inverse 21 ' 22 ' 2630 ; here, we choose the Moore-Penrose pseudoinverse 
which is well-defined for non-invertible matrices A and which we 
denote by A^ 1 . A brief review of the Moore-Penrose pseudoinverse 
is provided in the Supplementary Note. Applying the pseudoinverse 
to equation (5), we obtain 

^ = {k-^ 1 {-u^+^)- (9) 

Details of the derivation of equation (9) are given below in the 
Methods section. We emphasize that this pseudoinverse does not 
guarantee that (w^\u^)=0, explaining our previous remark on 
keeping this term intact. 

The steady-state density matrix p s = u 0 defined by Lp s = 0 is of 
particular interest. As a special case of equations (8) and (9), we can 

simplify the corrections k^ and to 

4 ;) =o, (10) 

P y=-ir 1 hp<!- 1) , (ii) 

see Methods for details. Corrections to k 0 and p s were previously 
derived 30 without using the Moore-Penrose pseudoinverse. The 
result for the density- matrix corrections in ref. 30 differs from ours 
[equation (11)] merely by a shift, 

pUUpW + ^M (12) 

where Cj is a constant. 

This shift can be interpreted as follows. Note that equation (5) does 
not have a unique solution since (L 0 — k^) is non-invertible. Indeed, 

for any given solution p^ of equation (5) its shifted counterpart from 
equation (12) is also a solution. The choice of a particular generalized 
inverse effectively selects a specific set of shift parameters Cj. The 
difference between the result in ref. 30 and ours merely reflects the 
different choices of generalized inverses. Since shifts of the form of 
equation (12) only affect the overall normalization of p s , our result 
for the steady state corrections is equivalent to that given in ref. 30. 
We properly normalize our result by defining 

P. = N[5^ 0 rfp«)] (13) 
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with N[A] =A/Tr [A] effecting normalization for any non-traceless 
matrix A. 

Finite-order PT truncates the series in equation (13) just as in 
closed-system PT. Let us denote the approximate result up to M-th 

order as p^ M = N \Y^ =0 oi j p { s j) ~\ . We next check whether p^ M 

indeed represents a proper density matrix, i.e., whether it is normal- 
ized, Hermitian and positive-semidefinite 15 . By virtue of the 
finite- order result p^' ,M is explicitly normalized. Hermiticity can be 
verified by noting that l^ 1 and fl_i map Hermitian operators to 
Hermitian operators. However, omission of higher-order terms in 
the truncation can render p a s ,M non-positive. In the examples pre- 
sented later on, this issue indeed occurs for certain parameter 
choices. Beyond mere conceptual concerns, non-positivity also pre- 
vents direct calculation of quantities such as state fidelity and 
von-Neumann entropy. This marks a key difference between 
closed-system PT and density-matrix PT. In closed-system PT, the 
approximate result is always a proper quantum state. By contrast, for 
density-matrix PT the approximate result may, strictly speaking, not 
be a proper density matrix. 

Similar issues with non-positivity of approximate density matrices 
are also encountered in quantum tomography, there caused by 
measurement errors. In this case, a maximum-likelihood method is 
typically used to reconstruct a physical density matrix from the non- 
positive approximation 33,34 . However, this method is impractical to 
apply to the large density matrices we are interested in. We next 
discuss an alternative strategy which reinstates positivity for large 
density matrices obtained within PT. 

Amplitude-matrix PT. We propose an amplitude-matrix PT to 
perturbatively construct an approximate steady-state density ma- 
trix which is manifestly positive. For this, recall that any Hermitian 
and positive-semidefinite matrix p can be decomposed in the 
form 35 : p = £C f - Following ref. 36, we will refer to f as the 
amplitude matrix. The decomposition p = ££ f is not unique: 
there are many choices for £ leading to the same matrix p. To 
eliminate these extra degrees of freedom, we choose £ to be lower 
triangular with real, non- negative diagonal elements. Existence and 
uniqueness of £ are then guaranteed by the Cholesky decomposi- 
tion. (In the case that p possesses a zero eigenvalue, the Cholesky 
decomposition is not unique but this issue can be bypassed by 
known strategies 37 ' 38 .) 

We start from the power series expansion of the steady- state 
amplitude matrix in a: £ s = y\_ n a j C^. Here, all matrices are 
again of lower triangular shape. Plugging this expansion into 
P s = CsCl and collecting terms of the same order in a, we obtain 

= (14) 



k=l 

Here, the super-operator Z 0 is defined via Z 0 A = Cj^A 1 " + A (cf^ 

The zero- order amplitude matrix ^ is directly obtained from equa- 
tion (14) by Cholesky decomposition and £^ is determined recur- 
sively from the system of linear equations in equation (15). We 
determine p(p in equation (15) by density- matrix PT; in this sense, 
the amplitude-matrix PT is based on the density- matrix PT. 

Once we truncate the amplitude matrix to M-th order, 
£a;M ^ ^ a^s^, we can compute the steady-state density matrix 



p^=^[C ;M (C ;M ) t ]- (16) 

Now, Ps am * s manifestly a proper density matrix: it is normalized, 
Hermitian and positive-semidefinite by construction. We note 
already that expectation values for observables obtained from ampli- 
tude-matrix PT, however, are not necessarily more accurate than 
those obtained from density- matrix PT, even if p a]M is slightly 
non -positive. The examples in the following section illustrate that 
the respective accuracies of density-matrix versus amplitude-matrix 
PT generally depend on the specific perturbation and system 
parameters. 

Amplitude-matrix PT is more involved if one or more eigenvalues 
of pf^ vanish, e.g. when p^ represents a pure state. In that case, Z 0 in 
equation (15) is non-invertible (see Methods) and thus a unique 
solution for equation (15) does not exist. Depending on the specific 
case, there may be infinitely many solutions or no solution. In the 
case of infinitely many solutions, we add a small identity- matrix 
component to p[° ) , i.e., p^ -> p^ + c 1 with 0 < c< 1 . The identity 
matrix then acts as a correction matrix which stabilizes the procedure 
of solving the linear equation [equation (15)] to provide a unique £^ . 
We utilize this correction matrix strategy in the second example 
discussed below. Whenever equation (15) has no solution, other 
forms of series expansions would need to be applied. We will not 
further consider that case in the present paper. The correction - 
matrix method and the validity of the series expansion are further 
discussed in the Methods section. 

Comparing PT with exact results. To illustrate the use and assess the 
accuracy of density-matrix and amplitude-matrix PT, we study two 
example systems, see Fig. 1. These two examples are neither claimed 
to be original nor of particular intrinsic interest. Our selection is 
guided by the intention to discuss two systems that are non-trivial, 
yet sufficiently small to still allow for a numerically exact treatment of 
the master equation. This way, we can compare steady- state 
expectation values from second-order PT with those obtained from 
exact diagonalization of L. We note again that the steady- state result 
obtained from finite-order density-matrix PT can be non-positive. 
This is indeed the case for some choices of parameters in the two 
following examples. Thus, we also apply the amplitude-matrix PT 
and compare results from the two perturbative treatments. 



a b 




Figure 1 | Schematic of the examples, (a) The system consists of a ring of 
four coupled spins with nearest-neighbor flip-flop interaction of strength t 
(dashed red lines). The spins are jointly driven with a coherent tone of 
strength e and frequency (blue arrow). Each spin is further subjected to 
local dissipation with spin relaxation rate y (green curly arrows), (b) The 
second system is composed of a three-resonator ring (black rectangles) 
coupled to a single qubit with Jaynes-Cummings interaction strength g 
(blue dotted line). The resonators are coupled to each other through 
photon hopping with rate k (red dashed lines). One of the resonators is 
driven coherently (blue arrow). Resonators and qubit are subject to local 
dissipation with photon decay rate y a and qubit relaxation rate y q (green 
curly arrows). 
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Unperturbed 




Kerr) 




Figure 2 | Results for four spins coupled in a ring. Shown are: the unperturbed result with respect to L 0 (blue dotted line), the exact result with respect to 
L (black solid line), the second-order density-matrix-PT result (red dashed line) and the second-order amplitude-matrix-PT result (green thin line). 

(a) and (b): | (g± ) | and (<7+ a± ) are plotted as a function of Sw/y for e/y = 0.8 and t/y = 0.4. The exact result is well approximated by the second-order 
PT result. The amplitude-matrix PT is slightly less accurate in this particular case which is more visible in (a* ). (c) and (d): The same is plotted except 
for t/y = 0.8. Although the deviation is relatively large, the shape of the exact result is still qualitatively captured by the second-order PT results. 



Ring of four coupled spins, driven and damped. We consider four 
spins coupled in a ring as shown in Fig. 1(a). The spins are coupled by 
flip-flop interaction with spin-spin coupling strength t. We assume 
all spins are driven equally with strength e and drive frequency co d . 
Within the Rotating Wave Approximation (RWA), the Hamiltonian 
describing this system is given by 

H= J2 [Sco o+a- +« +<)] 

(17) 

+ t E( <7 >«"+i+ hx -)- 

n 

Here, are the raising or lowering operators for the spin at site n, 
and dco = (o 0 — co d is the detuning between the spin frequency co 0 
and the drive frequency w d . Note that in equation (17), the time 
dependence of the drive has already been eliminated by working in 
the rotating frame. All four spins are coupled to a zero-temperature 
bath, leading to pure spin relaxation with relaxation rate y. Thus, the 
Liouville super-operator L is given by 

Lp=-i[H,p]+y^D[a-] p, (18) 

n 

where O [o~ ] p = a~ po+ - i <7+ a~ p - ^pa+ a~ is the usual 

dissipator for spin relaxation. 

As the perturbation, we now choose the spin-spin coupling and, 
thus, separate L into two parts, L = L 0 + 1 fl_i where 

l 0P = E { - 1 l SW < ^ + * {< + ^ ) A + 7 0 K" ] P } 

n 

describes the "atomic limit" in which the spin-spin coupling is 
absent, and the perturbation 

thp = tJ2-i[{<K + i+^°n + i)>P] (19) 

n 

captures the spin- spin coupling. Orders of perturbation theory are 
counted with respect to t. Even though the system possesses a high 



degree of symmetry, the steady states of the unperturbed and per- 
turbed systems are unique. Hence, the zero -eigenvalues of L 0 and fl- 
are non- degenerate and our non-degenerate PT for the steady state is 
applicable. 

We can thus implement second-order PT to compute the 
steady-state expectation values for several operators. We choose 
(jj - and the excitation number cj 1 + cr" since they fully capture 
the reduced density matrix of a single spin (which is the same 
for all spins due to symmetry). In Fig. 2, we compare results from 
density-matrix PT and amplitude-matrix PT to the exact and the 
unperturbed results. 

We first consider the case shown in Fig. 2(a) and (b) where the 
coupling strength t represents the smallest energy scale. Consistent 
with findings in ref. 4, we observe two symmetric resonance peaks in 
| (fff ) | for the unperturbed result (separate spins). When the coup- 
ling is present, the two peaks are shifted in position and become 
asymmetric, in agreement with the results from ref. 39. For 
( a i °\ )' we 0Dser ve a similar shift and asymmetry in the resonance 
peak due to the spin-spin coupling. Second-order PT nicely captures 
the above features. Note that the amplitude-matrix-PT result is actu- 
ally slightly less accurate than the density-matrix-PT result in this 
particular case. 

To illustrate the expected limitations of PT, we next increase the 
coupling strength so much that it matches the drive strength. Results 
for this parameter choice are shown in Fig. 2(c) and (d). 
Qualitatively, the shape of the curves from the exact calculation is 
still captured by the perturbative results. However, the results from 
PT show significant deviations from the exact result. Large deviations 
like this are not surprising since the "perturbation" parameter t now 
matches both e and y and PT breaks down. 

It is an interesting question whether a dimensionless parameter a 
can be identified that governs the convergence of this PT. To invest- 
igate this question, we recall that in simple cases of closed- system PT, 
a may be inferred from the expression for second-order corrections, 
a is then given by the ratio of the perturbation strength t and the 
difference between the relevant unperturbed eigenenergies; for 
instance, for the ground state we have: 
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Figure 3 | Eigenvalue spectrum of a driven-damped spin. The diagram 

shows the eigenvalue spectrum for system parameters chosen as 

dco/y = 0.2 and varying drive strength from e/y = 0 (red circles) to 
e/y = 0.25 (purple squares). The fourth eigenvalue (black diamond) 

corresponds to the steady state and is always zero. As illustrated by the blue 
trajectories, the flow of the eigenvalues is non-trivial and the eigenvalue(s) 
closest to zero switches from the complex pair to the purely real eigenvalue 
as e increases. 



t 



(closed system) 



(20) 



c 0 



^1 



For open systems, the unperturbed eigenvalues of the Hamiltonian 
are replaced by those of the super- operator L 0 . For the steady state, 
the relevant eigenvalue difference is that between the steady- state 
eigenvalue zero and that of the closest non-zero eigenvalue(s), 



(21) 



Even for a simple system like a driven-damped spin, the spectrum of 
complex eigenvalues depends on the system parameters 3co, e 
and y in a rather non- trivial way, see Fig. 3. The tuning of system 
parameters can even change the identity of the eigenvalue that is 
closest to zero. Hence, it is in general difficult to write the small 
parameter a in a simple form showing the dependence on the various 
system parameters. 

Qubit coupled to a driven and damped resonator ring. For our 

second example, we choose a system composed of a single qubit 
coupled to a three-resonator ring, see Fig. 1(b). The resonators are 
coupled among each other by photon hopping with rate k. The first 
resonator is driven with strength e and drive frequency co^. The qubit 
is coupled to the second resonator with a coupling strength g. Within 
RWA and in the frame co- rotating with the drive, the system 
Hamiltonian H is given by 



H= ^2 fa a\a n + K(a\a n+l + fl„4+i)] +e(ai+a\ 

n 

+ 3co a + g~ +g(a 2 G + +a\o~^j . 



(22) 



Here, a n (aj) is the annihilation (creation) operator for photons in 
the resonator at site n and 3w = w — coj is the detuning between the 
bare resonator frequency co and the drive frequency co^. The qubit is 
assumed in resonance with the bare resonator frequency. Qubit and 
resonators are each coupled to a zero-temperature bath, inducing 



qubit relaxation and photon decay with rates y q and y a , respectively. 
The super- operator L is thus given by 



Lp- 



i[H,p}+y a J2®[an]p + Y q B[°-}p- 



(23) 



The qubit- resonator coupling mediates two effects: an indirect co- 
herent drive on the qubit, and correlations between the resonator 
ring and the qubit. We wish to treat the correlation effects 
perturbatively. To do so, we separate the two effects by means of a 
coherent displacement as follows. Note that in the absence of qubit - 
resonator coupling, the drive places the eigenmodes /n= 1, 2, 3 of the 

resonator ring in coherent states with amplitudes (a^)^ 



oc^ where 



6 fx ,1 2n ^ 
oco + Ik cos 

V3 V 3 



expl i 



2nfi 



(24) 



and 0^ = a n exp ^ Lir^J is the annihilation operator for 

photons in mode fi. We thus displace according to 



(25) 



and rewrite the Liouville super- operator as 

lp=-i[H 0 +gH 1 ,p}+y a ^B[a' >l ]p+y q BlcT-}p. 



Here, H 0 is the unperturbed Hamiltonian 

H 0 = ( 3co + 2k cos — ^ 



(26) 



+ 3tOO + G + (e e ff<7 + +6*ff°" 



in which the displaced eigenmodes a'^ and the qubit remain 
decoupled. One finds that the effective drive on the qubit is given 

by e e ff = -7= 5^ exp ( — i — ^ ) . The perturbation gH 1 describes 

the remaining coupling between the displaced eigenmodes and the 
qubit, 



3 / 



(27) 



The perturbative treatment at the level of the master equation simply 
follows this separation and decomposes L = L 0 +^Li into the 
unperturbed super-operator 



L 0 p = 



! [Ho,p]+y fl $3D[5J,]p+y,D[ff-]p, 



(28) 



and a perturbation which only captures the remaining coupling, 

ghp=-i\gH l9 p\. (29) 

The order of perturbation theory here is counted with respect to g. 

The steady state of D_ 0 is a product state composed of a density 
matrix for the resonator ring and one for the qubit. The resonator 
ring is in a pure state (all displaced eigenmodes in the vacuum state). 
As a consequence, the unperturbed density matrix has multiple 
eigenvalues zero. Therefore, when implementing amplitude-matrix 
PT, we employ the correction matrix method mentioned above (with 
parameter c= 10" 9 ). 

We will focus on expectation values at site 1, specifically (ai) and 
(tti ) = (a\ai ), as a function of the drive frequency, expressed in terms 
of the detuning 3w. In the unperturbed, i.e., decoupled case (g = 0), 
we expect two resonances at 3co = —2k and 3co = k corresponding 
to the eigenmodes of the resonator ring. Once the qubit is coupled to 
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Figure 4 | Results for single qubit coupled to a resonator ring. The color scheme of Fig. 2 is used. To apply the amplitude-matrix PT, the correction 
matrix is employed with the parameter c = 10" 9 . We plot l(a x ) I , and as a function of Sw/e for k/c = 10, g/e = 0.5 and y q /^ = y a / e = 0.05. The 

resonance at Sco = 0 is well captured by second-order PT. The amplitude-matrix PT and density-matrix PT perform with nearly identical accuracy. 



the resonator ring, we expect a resonance at dco = 0 originating from 
the qubit's response to the drive. We monitor this response by cal- 
culating the expectation value of o~. 

As shown in Fig. 4, we confirm that the resonance at dco = 0, a key 
consequence of the coupling, is successfully captured by second- 
order PT. Specifically, we consider the case of drive and coupling 
strengths e and g small compared to the hopping rate k but large 
compared to the relaxation and decay rates y q and y a . Note that 
the perturbation parameter g is not the smallest energy scale in this 
case; nonetheless, PT holds. The expectation values of a l5 n x and o~ 
are shown in Fig. 4(a), (b) and (c) respectively. The results from 
second-order PT are in good agreement with the exact result. 
Amplitude-matrix PT and density- matrix PT here give results with 
nearly identical accuracy. The saturation effect visible in Fig. 4(c) 
shows that the qubit is in the nonlinear regime. These results also 
illustrate that the correction matrix method required for the ampli- 
tude-matrix PT is reliable and yields results consistent with the exact 
solution. 



Discussion 

We have detailed a perturbative approach to Markovian open 
quantum systems and developed a non- degenerate PT based on 
the Lindblad master equation. This density- matrix PT recursively 
determines corrections to the eigenvalues and right eigenstates of 
the Liouville super-operator L. As a result of finite-order truncation, 
such a perturbative scheme may lead to non-positive steady-state 
"density matrices" which prevent direct calculation of quantities like 
the state fidelity. The issue of non-physical states, which does not 
occur in closed-system PT, can be tackled by a modified perturbative 
scheme based on the amplitude matrix. With two example systems, 
we have illustrated that the approximate results are in excellent 
agreement with exact results for representative parameter choices. 
The expectation values obtained from density- matrix PT showed 
good agreement in the two examples even if the truncated density 
matrix was slightly non-positive. 

The perturbative treatment presented here is suitable for systems 
of sizes that cannot be handled by exact solution of the quantum 
master equation. An interesting future application of this PT thus 
consists of the study of open quantum systems with lattice structure, 
such as the open Jaynes-Cummings lattice. Promising experimental 
progress 5 indicates that such open lattices can indeed be realized in 
the circuit QED architecture, and may serve as valuable open-system 
quantum simulators 11 . Openness and relatively large size of such 
systems make the theoretical investigation challenging. We believe 
that the developed open- system PT will provide a useful tool in 
exploring the physics of open lattice systems. 

Methods 

Perturbative corrections. We wish to prove that the expression of in equation 
(9), rewritten as 



««=(lo-4 0) ) a. 

is indeed a solution to equation (5), i.e., it satisfies 



(30) 



(31) 



where the operators^ are defined as^ = — D_ i m jj 1 ^ + . Solving equation ( 3 1 ) for 
uf is a standard linear algebra problem. The necessity for working with a 
pseudoinverse lies in the fact that (Lq — X^) is singular and non-Hermitian. This 
prevents us from using the ordinary inverse to solve for ufjfi . Here, we employ the 
Moore-Penrose pseudoinverse denoted by "~ 1 . ( A brief review of the pseudoinverse is 
given in the accompanying Supplementary Note.) While this choice is not unique, it is 
convenient since the Moore-Penrose pseudoinverse can be computed efficiently via 
singular value decomposition. 

After plugging i/jp from equation (30) into equation (31), it is clear that the proof 
amounts to verifying that 



(l 0 -i(°>)(l„-i(»>) 



(32) 



Employing the defining properties of the Moore-Penrose pseudoinverse, one finds 
that the claim of equation (32) can be written in the equivalent form 



(33) 



where P L _ ^o) is the orthogonal projector onto the range of (L 0 — Since every 
projector acts as the identity on vectors from its range, equation (33) holds iffj M 
belongs to the range 9ft At of D =D Lo _ /L (o) . Furthermore, since ^ Lo _^°) i s an orthogonal 
projector, the range 3^ and the null space 91^ of are orthogonal subspaces. 

This implies that^ belongs to 9^ if and only if is orthogonal to 91^. To see that 
fjn-Lyi^, note that (w^ ,fj^) = 0 as a direct consequence of equation (5). To prove that 
f^X-tyl^ we now simply make use of the fact that is spanned by the left eigenstate, 
i.e., 9T At = span{w/ 0) }. 

We verify the last statement as follows. Since is the left eigenstate of L 0 with 
eigenvalue ^ 0) , it is clear that (wf\(l 0 - A^)A) = 0 for any matrix A. Thus, is 
orthogonal to the range of (L 0 — Since (L 0 — X^) and PY ; (°) share the same 
range, is also orthogonal to the range 9^ of P 1q _ ; ( 0 ) . Thus, is an element of 
the null space 91^ of _ ; ( 0 ) . Assuming that X^ is a non-degenerate eigenvalue of D_ 0 , 
we thus find that 9T M is spanned by , as stated. 

Steady-state corrections. For the steady state p s defined by Lp s = 0, i.e. p s = u 0 , the 
recursion relations (8) and (9) can be simplified significantly. To see this, recall that L 
and L 0 are proper generators of the quantum dynamical semigroup. In order to be 
trace preserving (a necessary condition for a proper generator), the identity 1 must be 
the left eigenstate of L and L 0 with eigenvalue zero, i.e. w 0 = w ® = 1. It follows that 1 
is also the left eigenstate of D_i with eigenvalue zero since Li = L — L 0 and thus 
TrfO-ix] = 0 for any operator x. It is straightforward to show from Tr[lLix] = 0 and 
equations (8) and (9) that VjeN, 



(34) 
(35) 



Series expansion of amplitude matrices. The series expansion of £ s is valid if all £^ 
can be determined according to 
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ZcCM-ge^y, (36) 
k=l 

which is equation (15). Here Z 0 was defined via Z 0 A = ^ 0) A t +A(Cf ) ) t and } is 
determined through Cholesky decomposition: =Cs°^(Cs°^- Equation (36) is a 
system of linear equations and thus has a unique solution if and only if Z 0 is invertible. 

Whether Zq is invertible depends on the form of Q° > which in turn depends on 
the Hermitian and positive -semidefinite p[°\ If one of the eigenvalues of p^ is 
zero, there is a corresponding eigenvector \\J/) such that p^\\J/) =0. Consider the 
decomposition: p^ = hh where h is a Hermitian matrix. The existence of this 
decomposition can be simply proven by writing p^ in its eigen-decomposition 
form. Since is the eigenvector of p^ with eigenvalue zero, is also the 
eigenvector of h with eigenvalue zero, i.e. h = 0. Moreover, due to the fact that 
C s (0) (C s (0) ) t =pf ] =hh y and h are unitarily right equivalent 40 , i.e. ^ 0) =hS where 
S is a unitary matrix. Thus, \\Jj) is also the left eigenvector of £^ with eigenvalue 
zero, i.e. (£^)^|iA) =S ,f h\iJ/) =0. Now, there must be a right eigenvector of C^K 
denoted by |0), that corresponds to the same eigenvalue (which is zero), i.e. 
Cs°^|0) =0- We can show by directly substitution that Z o (T|0)(0|) = 0 where T is 
the matrix corresponding to Gaussian elimination which transforms (|0) (0|) to a 
lower triangular matrix. Therefore, Zq is not invertible if p^ contains at least one 
eigenvalue zero. 

If Zq is not invertible, equation (36) either has infinitely many solutions or no 
solution. The former case, in which there are infinitely many solutions, can be 
bypassed if we shift the eigenvalues of p^ away from zero. To do so, we consider a 
shift by an identity- matrix component according to 

pW^pW+d. (37) 

We choose the auxiliary parameter c as close to zero as possible while maintaining 
numerical stability. In this way, we obtain a procedure to obtain a unique solution to 
equation (36). This method is similar to the correction matrix method 37 ' 38 for the 
Cholesky decomposition of matrices with eigenvalue(s) zero. There, a small 
diagonal correction matrix is also added to the original matrix to avoid the 
eigenvalue(s) zero. 

If we encounter the case in which there is no solution, we cannot determine . 
In fact, a two-level system coupled to finite temperature bath with D [a + ] term treated 
as the perturbation belongs to this case. The failure to determine £^ indicates that the 
series expansion of [ s is invalid. This originates from the fact that if p^ contains 
any zero eigenvalues, the leading order term of some elements of £ s may be of the 
order of a 112 (or a 3/2 , etc.) instead of a 0 . For real functions, an analogy would be the case 
y(a) = x 2 (a) where;/ can be written as a power series in a. If the leading order term of v 
is of the order of a, x is a series that only contains half- integer orders of a and thus is 
not a proper Taylor series. Different expansion types would be needed in this case 
which we do not discuss further in the present paper. 
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